---
title: "Analysis for: The Dynamics of Electoral Manipulation and Institutional Trust in Democracies (POQ)"
author: "M. Higashijima, H. Kadoya, and Y. Yanai"
date: "2023-11-12"
output: pdf_document
---


```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r global_option, include = FALSE}
knitr::opts_chunk$set(
  echo = TRUE,      
  warning = FALSE,  
  message = FALSE,  
  fig.width = 6,    
  fig.height = 3)   
if (knitr::is_latex_output()) {
  knitr::opts_chunk$set(dev = "cairo_pdf")
}
```

```{r, include = FALSE}
pacman::p_load(tidyverse, 
               estimatr, 
               MASS,
               interplot,
               patchwork,
               texreg,
               data.table,
               ggrepel,
               sjPlot,
               interactions,
               tinytex)
```

Function to show marginal effects of interaction models.
```{r}
meplot <- function(model,
                   var1,
                   var2,
                   data,
                   int,
                   vcov,
                   ci = 0.95,
                   xlab = var2,
                   ylab = paste("Marginal Effect of", var1)) {
  alpha <- 1 - ci
  z <- qnorm(1 - alpha / 2)
  beta_hat <- coef(model)
  cov <- vcov
  x0 <- pull(data, var2)
  z0 <- seq(min(x0, na.rm = TRUE), 
            max(x0, na.rm = TRUE),
            length.out = 1000)
  dy_dx <- beta_hat[var1] + beta_hat[int] * z0
  se_dy_dx <- sqrt(cov[var1, var1] + z0^2 * cov[nrow(cov), ncol(cov)] +
                     2 * z0 * cov[var1, ncol(cov)])
  upr <- dy_dx + z * se_dy_dx
  lwr <- dy_dx - z * se_dy_dx
  ggplot(data = NULL, aes(x = z0, y = dy_dx)) +
    labs(x = xlab, y = ylab) +
    geom_hline(yintercept = 0,
               linetype = "dashed",
               color = "black") +
    geom_ribbon(aes(ymin = lwr, ymax = upr),
                alpha = 0.3) +
    geom_line(aes(z0, dy_dx),
              color = "black")
}
```


# Data

## Data prep

Read individual-level data.
```{r}
IND_nd <- read_rds("data_individual_HKY_POQ_2023.Rds")
```

Subset of the data: democracy
```{r}
IND_nd_dem <- IND_nd |> 
  filter(democracy == 1)
```

Subset of the data: complete cases.
```{r}
IND_nd_dem_sub <- IND_nd_dem |> 
  filter(!is.na(early),
         !is.na(dif_days10),
         !is.na(fraud2),
         !is.na(nelda17),
         !is.na(nelda18),
         !is.na(gdpgrowth_new),
         !is.na(gdppc_change),
         !is.na(margin_pty_lag),
         !is.na(age),
         !is.na(sex),
         !is.na(employment),
         !is.na(urban))
```


Clean up some variables.
```{r}
MYD <- IND_nd_dem |> 
  mutate(female = ifelse(sex == 1, 1, 0),
         employment2 = ifelse(survey_wave == "Afro1",
                              employment + 1,
                              employment),
         unemployment = ifelse(employment2 == 2, 1, 0)) |> 
  mutate(nelda17 = ifelse(nelda17 == "yes", 1, 0))
```

Data used in regression analyses.
```{r}
MYDout <- MYD |> 
  dplyr::select(cowcode, year, survey_wave,
                GOV, PMT, ETN,
                early, dif_days10, nelda17, nelda18, fraud,
                fraud1, fraud2,
                gdpgrowth_new, margin_pty_lag,
                age, female, unemployment, urban,
                pid_ruling,
                pid_ruling_voting,
                coalition, win, gov1me, gov2me, gov3me,
                v2x_polyarchy_excfraud, v2x_libdem_excfraud,
                date_dif,
                gdppc_change
                ) |> 
  mutate(GOV = GOV + 5,
         PMT = PMT + 5,
         ETN = ETN + 5,
         fraud = fraud + 1)
```


Subset of the data: complete cases.
```{r}
MYDout2 <- MYDout |> 
  filter(!is.na(early),
         !is.na(dif_days10),
         !is.na(fraud2),
         !is.na(nelda17),
         !is.na(nelda18),
         !is.na(gdpgrowth_new),
         !is.na(gdppc_change),
         !is.na(margin_pty_lag),
         !is.na(age),
         !is.na(female),
         !is.na(unemployment),
         !is.na(urban)) |> 
  mutate(cowcode = factor(cowcode),
         year = factor(year))
#MYDout2 <- read_rds("data_reg_HKY_POQ_2023.Rds")
```


## Descriptive stats

Figure B2: Distribution of blatant electoral fraud.
```{r}
hist_fraud <- IND_nd |> 
  dplyr::select(cowcode, year, fraud2, democracy) |> 
  distinct() |> 
  filter(!is.na(democracy)) |> 
  mutate(democracy = factor(democracy,
                            levels = c(0, 1),
                            labels = c("Autocracy", "Democracy"))) |> 
  ggplot(aes(x = fraud2)) +
  geom_histogram(color = "black",
                 boundary = 1) +
  labs(x = "Blatant fraud") +
  facet_grid(row = vars(democracy))
plot(hist_fraud)
```


Figure B3: Scatterplot of blatant fraud vs polyarchy among democracies.
```{r}
D_fig <- IND_nd_dem_sub |> 
  filter(!is.na(GOV)) |> 
  dplyr::select(cowcode, country, year, fraud2, v2x_polyarchy_excfraud) |> 
  distinct() |> 
  mutate(name = paste(country, year, sep = "-"))
fraud_pol_dem_lab <- 
  ggplot(D_fig, aes(x = v2x_polyarchy_excfraud, y = fraud2)) +
  geom_point(color = "gray40") +
  labs(x = "Polyarchy", y = "Blatant fraud") +
  xlim(0, 1) + ylim(0, 1) +
  geom_text_repel(data = subset(D_fig, fraud2 > 0.65),
                  aes(v2x_polyarchy_excfraud, fraud2, label = name),
                  hjust = 1.2,
                  color = "black",
                  segment.color = "gray")
plot(fraud_pol_dem_lab)
```

Figure 1: Popular trust in the government, parliament, and EMBs.
```{r}
sum(!is.na(MYDout2$GOV))
dist_gov <- MYDout2 |> 
  filter(!is.na(GOV)) |> 
  mutate(Gov = factor(GOV,
                      labels = c("1. Not at all",
                                 "2. Just a little",
                                 "3. Somewhat",
                                 "4. A lot"))) |> 
  ggplot(aes(x = Gov, y = (..prop..) * 100, group = 1)) +
  geom_bar(color = "black",
           stat = "count") +
  labs(x = "Trust in government\n(N = 335,453)", 
       y = "percent") + 
  ylim(0, 40) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
#plot(dist_gov)

sum(!is.na(MYDout2$PMT))
dist_pmt <- MYDout2 |> 
  filter(!is.na(PMT)) |> 
  mutate(Pmt = factor(PMT,
                      labels = c("1. Not at all",
                                 "2. Just a little",
                                 "3. Somewhat",
                                 "4. A lot"))) |> 
  ggplot(aes(x = Pmt, y = (..prop..) * 100, group = 1)) +
  geom_bar(color = "black",
           stat = "count") +
  labs(x = "Trust in parliament\n(N = 325,289)",
       y = "percent") +
  ylim(0, 40) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
#plot(dist_pmt)

sum(!is.na(MYDout2$ETN))
dist_emb <- MYDout2 |> 
  filter(!is.na(ETN)) |> 
  mutate(Etn = factor(ETN,
                      labels = c("1. Not at all",
                                 "2. Just a little",
                                 "3. Somewhat",
                                 "4. A lot"))) |> 
  ggplot(aes(x = Etn, y = (..prop..) * 100, group = 1)) +
  geom_bar(color = "black",
           stat = "count") +
  labs(x = "Trust in EMBs\n(N = 105,382)",
       y = "percent") +
  ylim(0, 40) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
#plot(dist_emb)
plot(dist_gov | dist_pmt | dist_emb)
```

Days since the last election.
```{r}
DDD <- IND_nd_dem_sub |> 
  dplyr::select(survey_wave, country, dif_days) |> 
  distinct() 

mean(DDD$dif_days)
sd(DDD$dif_days)
```

Figure 2: The number of days elapsed since the last legislative election before the survey.
```{r}
DD <- DDD |> 
  ggplot(aes(dif_days)) +
  geom_histogram(color = "black",
                 boundary = 0) +
  labs(x = "Days since the last election")
DD
```


List of country-survey-year (election year) included in the dataset
```{r}
list_elections <- IND_nd_dem |> 
  dplyr::select(country, 
                survey_wave,
                year,
                early,
                fraud2) |> 
  distinct() |> 
  arrange(country, year)
list_elections
```


Descriptive Stat Table
```{r}
Individual <- data.table(
  Gov = MYDout2$GOV,
  Pmt = MYDout2$PMT,
  EMB = MYDout2$ETN,
  Age = MYDout2$age,
  Female = MYDout2$female,
  Unemployment = MYDout2$unemployment,
  Urban = MYDout2$urban
)
```

Country-election data.
```{r}
MYD2 <- MYDout2 |> 
  filter(!is.na(GOV))
CED <- MYDout2 |>
  filter(!is.na(GOV)) |> 
  dplyr::select(cowcode, year, #survey_wave,
                date_dif,
                early, dif_days10, fraud2, 
                v2x_polyarchy_excfraud, v2x_libdem_excfraud,
                nelda17, nelda18,
                gdpgrowth_new, gdppc_change, margin_pty_lag, coalition) |> 
  distinct() |> 
  mutate(CE_ID = 1:n())

CYlevel <- data.table(
  Early = CED$early,
  Days = CED$dif_days10,
  Fraud = CED$fraud2,
  Polyarchy = CED$v2x_polyarchy_excfraud,
  `Liberal Democracy` = CED$v2x_libdem_excfraud,
  NELDA17 = CED$nelda17,
  NELDA18 = CED$nelda18,
  GDP = CED$gdpgrowth_new, 
  `GDP (per capita) change` = CED$gdppc_change,
  Margin = CED$margin_pty_lag,
  Coalition = CED$coalition
)
```

LaTeX table of Stats (Table B1)
```{r}
datasets <- list("Individual-level" = Individual, 
                 "Election-level" = CYlevel)
variables <- list(names(Individual), names(CYlevel))
labels <- list(names(Individual), names(CYlevel))
colnames <- c("Mean", "Median", "SD", "Min", "Max", "N")

# descriptive function
my_descriptives = function(x) {
    x = as.numeric(x)
    m = mean(x, na.rm = TRUE)
    md = median(x, na.rm = TRUE)
    sd = sd(x, na.rm = TRUE)
    min = min(x, na.rm = TRUE)
    max = max(x, na.rm = TRUE)
    n = as.integer(sum(!is.na(x)))
    return(c(m, md, sd, min, max, n))
}
source("https://gist.githubusercontent.com/sdaza/c4d1089a501d3567be9fb784b1c5a6ab/raw/7a6237662d4ca2abcc4941ebbf21aa23d89721fc/create_descriptive_tables.R")
createDescriptiveTable(datasets,
    summary_function = my_descriptives,
    column_names = colnames,
    variable_names = variables,
    variable_labels = labels,
    arraystretch = 1.3,
    title = "Descriptive statistics",
    label = "tab:descriptive"
    #file = "figs/descriptives.tex" 
    )
```

```{r}
MYD |> 
  filter(!is.na(ETN)) |> 
  dplyr::select(ETN, early, dif_days10, nelda17, nelda18, CEN,
                gdpgrowth_new, margin_pty_lag,
                age, female, unemployment, urban) |> 
  summary()
```


Figure B1: Distribution of the timing of early elections relative to the scheduled date.
```{r}
hist_days_early <- CED |> 
  filter(early == 1) |> 
  dplyr::select(date_dif, year, cowcode) |> 
  distinct() |> 
  ggplot(aes(x = date_dif)) +
  geom_histogram(binwidth = 50, 
                 color = "black",
                 boundary = -0) +
  labs(x = "Days between the election and the term expiration")
hist_days_early
```

List of countries analyzed (Table B2)
```{r}
countries_inReg <- MYDout2 |>
  filter(!is.na(GOV)) |> 
  dplyr::select(cowcode) |> 
  mutate(cowcode = as.character(cowcode)) |> 
  distinct() 
countries_inReg <- IND_nd |> 
  dplyr::select(country, cowcode) |> 
  mutate(cowcode = as.character(cowcode)) |> 
  distinct() |> 
  right_join(countries_inReg, by = "cowcode") |> 
  arrange(country)
countries_inReg$country
```


# Regression: Effect of election timing

## Models

Models
```{r}
m_gov <- alist(
  m1 = GOV ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud + 
    age + female + unemployment + urban,
  m2 = GOV ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + 
    age + female + unemployment + urban,
  m3 = GOV ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud + 
    fraud2 + 
    age + female + unemployment + urban,
  m4 = GOV ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud + 
    fraud2 + 
    age + female + unemployment + urban,
  m5 = GOV ~ early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud + fraud2 + 
    age + female + unemployment + urban,
  m6 = GOV ~ early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + fraud2 + 
    age + female + unemployment + urban) |> 
  map(.f = formula)

m_gov_f <- alist(
  m1 = GOV ~ 0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud + 
    age + female + unemployment + urban + 
    cowcode + year,
  m2 = GOV ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + 
    age + female + unemployment + urban + 
    cowcode + year,
  m3 = GOV ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud + 
    fraud2 + 
    age + female + unemployment + urban + 
    cowcode + year,
  m4 = GOV ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud + 
    fraud2 + 
    age + female + unemployment + urban + 
    cowcode + year,
  m5 = GOV ~  0 + early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud +  fraud2 + 
    age + female + unemployment + urban + 
    cowcode + year,
  m6 = GOV ~  0 + early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + fraud2 + 
    age + female + unemployment + urban + 
    cowcode + year) |> 
  map(.f = formula)

m_pmt <- alist(
  m7 = PMT ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud + 
    age + female + unemployment + urban,
  m8 = PMT ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + 
    age + female + unemployment + urban,
  m9 = PMT ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud + 
    fraud2 + 
    age + female + unemployment + urban,
  m10 = PMT ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud + 
    fraud2 + 
    age + female + unemployment + urban,
  m11 = PMT ~ early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud +  fraud2 + 
    age + female + unemployment + urban,
  m12 = PMT ~ early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + fraud2 + 
    age + female + unemployment + urban) |> 
  map(.f = formula)

m_pmt_f <- alist(
  m7 = PMT ~ 0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud + 
    age + female + unemployment + urban + 
    cowcode + year,
  m8 = PMT ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + 
    age + female + unemployment + urban + 
    cowcode + year,
  m9 = PMT ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud + 
    fraud2 + 
    age + female + unemployment + urban + 
    cowcode + year,
  m10 = PMT ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud + 
    fraud2 + 
    age + female + unemployment + urban + 
    cowcode + year,
  m11 = PMT ~  0 + early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud + 
    age + female + unemployment + urban + 
    cowcode + year,
  m12 = PMT ~  0 + early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud +
    age + female + unemployment + urban + 
    cowcode + year) |> 
  map(.f = formula)

m_emb <- alist(
  m13 = ETN ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud + 
    age + female + unemployment + urban,
  m14 = ETN ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + 
    age + female + unemployment + urban,
  m15 = ETN ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud + 
    fraud2 + 
    age + female + unemployment + urban,
  m16 = ETN ~ early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud + 
    fraud2 + 
    age + female + unemployment + urban,
  m17 = ETN ~ early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud +  fraud2 + 
    age + female + unemployment + urban,
  m18 = ETN ~ early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + fraud2 + 
    age + female + unemployment + urban) |> 
  map(.f = formula)

m_emb_f <- alist(
  m13 = ETN ~ 0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud + 
    age + female + unemployment + urban + 
    year,
  m14 = ETN ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud + 
    age + female + unemployment + urban + 
    year,
  m15 = ETN ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud + 
    fraud2 + 
    age + female + unemployment + urban + 
    year,
  m16 = ETN ~  0 + early * dif_days10 + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud + 
    fraud2 + 
    age + female + unemployment + urban + 
    year,
  m17 = ETN ~  0 + early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_polyarchy_excfraud +  fraud2 + 
    age + female + unemployment + urban + 
    year,
  m18 = ETN ~  0 + early * dif_days10 * coalition + 
    nelda17 + nelda18 + gdpgrowth_new + gdppc_change + margin_pty_lag +
    v2x_libdem_excfraud +  fraud2 + 
    age + female + unemployment + urban + 
    year) |> 
  map(.f = formula)
```

Outcome: Government trust
```{r}
fits_gov_f <- m_gov_f |> 
  map(.f = lm,
      data = MYDout2)
fits_gov <- m_gov |> 
  map(.f = lm_robust, 
      data = MYDout2, 
      clusters = cowcode,
      fixed_effects = ~ cowcode + year,
      se_type = "stata")
res_gov <- fits_gov |> 
  map(.f = tidy) |> 
  bind_rows() |> 
  filter(term %in% c("early", "early:dif_days10")) |> 
  dplyr::select(outcome, term, estimate, p.value:conf.high) |> 
  print()
```

Outcome: Parliament trust
```{r}
fits_pmt_f <- m_pmt_f |> 
  map(.f = lm,
      data = MYDout2)
fits_pmt <- m_pmt |> 
  map(.f = lm_robust, 
      data = MYDout2, 
      clusters = cowcode,
      fixed_effects = ~ cowcode + year,
      se_type = "stata")
res_pmt <- fits_pmt |> 
  map(.f = tidy) |> 
  bind_rows() |> 
  filter(term %in% c("early", "early:dif_days10")) |> 
  dplyr::select(outcome, term, estimate, p.value:conf.high) |> 
  print()
```

Outcome: EMB trust
```{r}
fits_emb_f <- m_emb_f |> 
  map(.f = lm,
      data = MYDout2)
fits_emb <- m_emb |> 
  map(.f = lm_robust, 
      data = MYDout2, 
      clusters = cowcode,
      fixed_effects = ~ year,
      se_type = "stata")
res_emb <- fits_emb |> 
  map(.f = tidy) |> 
  bind_rows() |> 
  filter(term %in% c("early", "early:dif_days10")) |> 
  dplyr::select(outcome, term, estimate, p.value:conf.high) |> 
  print()
```


## Interaction effect

Without FRAUD variable.
Figure A1: The marginal effect of the early election on the popular trust of the government,
parliament, and EMBs, conditional on the number of days since the last election,
without controlling for blatant electoral fraud.
```{r}
intp_gov1 <- meplot(fits_gov_f$m1, 
   vcov = vcov(fits_gov$m1), 
   data = MYD,
   var1 = "early", 
   var2 = "dif_days10",
   int = "early:dif_days10") +
  labs(x = "Days since last election",
       y = "Effect of ealry election",
       title = "Government") +
  ylim(-1, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intp_gov1)

intp_pmt1 <- meplot(fits_pmt_f$m7, 
   vcov = vcov(fits_pmt$m7), 
   data = MYD,
   var1 = "early", 
   var2 = "dif_days10",
   int = "early:dif_days10") +
  labs(x = "Days since last election",
       #y = "Effect of ealry election",
       y = "",
       title = "Parliament") +
  ylim(-1, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intp_pmt1)

intp_emb1 <- meplot(fits_emb_f$m13, 
   vcov = vcov(fits_emb$m13), 
   data = MYD,
   var1 = "early", 
   var2 = "dif_days10",
   int = "early:dif_days10") +
  labs(x = "Days since last election",
       #y = "Effect of ealry election",
       y = "",
       title = "EMB") +
  ylim(-1, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intp_emb1)
plot(intp_gov1 | intp_pmt1 | intp_emb1)
```

With FRAUD variable.
Figure 3: The marginal effect of the early elections on popular trust in the government
(left panel), parliament (middle), and EMBs (right), conditional on the number of days
elapsed since the last election.
```{r}
intp_gov2 <- meplot(fits_gov_f$m3, 
   vcov = vcov(fits_gov$m3), 
   data = MYD,
   var1 = "early", 
   var2 = "dif_days10",
   int = "early:dif_days10") +
  labs(x = "Days since last election",
       #x = "",
       y = "Effect of ealry election",
       title = "Government") +
  ylim(-1, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intp_gv2)

intp_pmt2 <- meplot(fits_pmt_f$m9, 
   vcov = vcov(fits_pmt$m9), 
   data = MYD,
   var1 = "early", 
   var2 = "dif_days10",
   int = "early:dif_days10") +
  labs(x = "Days since last election",
       y = "",
       #y = "Effect of early election",
       title = "Parliament") +
  ylim(-1, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intp_pmt2)

intp_emb2 <- meplot(fits_emb_f$m15, 
   vcov = vcov(fits_emb$m15), 
   data = MYD,
   var1 = "early", 
   var2 = "dif_days10",
   int = "early:dif_days10") +
  labs(x = "Days since last election",
       #x = "",
       y = "",
       #y = "Effect of ealry election",
       title = "EMB") +
  ylim(-1, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intp_emb2)

plot(intp_gov2 + intp_pmt2 + intp_emb2)
```



## Predicted values: Three-way interaction

Government trust
```{r}
gov_coalition <- plot_model(fits_gov_f$m5, 
           type = "pred", 
           terms = c("dif_days10", "early", "coalition"),
           vcov.fun = vcov(fits_gov$m5)) +
  ylim(1, 3) +
  labs(x = "Days since last election",
       y = "Government",
       title = "Predicted values of trust") +
  scale_color_brewer(palette = "Set1",
                     name = "",
                     labels = c("As scheduled", "Early"),
                     guide = "none") +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
```


Parliament
```{r}
pmt_coalition <- plot_model(fits_pmt_f$m11, 
           type = "pred", 
           terms = c("dif_days10", "early", "coalition"),
           vcov.fun = vcov(fits_pmt$m11)) +
  ylim(1.5, 3.5) +
  labs(x = "Days since last election",
       y = "Parliament",
       title = "") +
  scale_color_brewer(palette = "Set1",
                     name = "",
                     labels = c("As scheduled", "Early")) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10) +
  theme(legend.position="bottom")
```

Figure A2: Predicted values of trust in government and parliament with the three-way
interaction.
```{r}
plot(gov_coalition / pmt_coalition)
```


## Table

Main results: Table 1
```{r, results = 'asis'}
texreg(list(fits_gov$m3,
            fits_pmt$m9,
            fits_emb$m15),
   stars = numeric(0),
   override.se = list( # p-values instead of se
     fits_gov$m3$p.value,
     fits_pmt$m9$p.value,
     fits_emb$m15$p.value
   ),
   include.ci = FALSE,
   include.rsquared = FALSE,
   include.adjrs = TRUE,
   include.nobs = TRUE,
   include.fstatistic = FALSE,
   include.rmse = FALSE,
   digits = 3,
   custom.model.names = c("gov", "pmt", "emb"),
   custom.note = "Note: p-values calculated with robust standard errors clustered by country are in parentheses. p-values refer to two-tailed tests. The first two models include year-fixed effects and all models include country-fixed effects.",
   caption = "",
   caption.above = TRUE,
   label = "tab:ols-res-main",
   float.pos = "t",
   #file = "figs/table_ols_res_main.tex",
   dcolumn = FALSE)
```

Government: Table A1
```{r, results = 'asis'}
texreg(list(fits_gov$m1,
            fits_gov$m2,
            fits_gov$m3,
            fits_gov$m4,
            fits_gov$m5,
            fits_gov$m6),
   stars = numeric(0),
   override.se = list( # p-values instead of se
     fits_gov$m1$p.value,
     fits_gov$m2$p.value,
     fits_gov$m3$p.value,
     fits_gov$m4$p.value,
     fits_gov$m5$p.value,
     fits_gov$m6$p.value),
   include.ci = FALSE,
   include.rsquared = FALSE,
   include.adjrs = TRUE,
   include.nobs = TRUE,
   include.fstatistic = FALSE,
   include.rmse = FALSE,
   digits = 3,
   custom.model.names = paste("Model", 1:6),
   custom.note = "Note: p-values calculated with robust standard errors clustered by country are in parentheses. p-values refer to two-tailed tests. All models include country and year fixed-effects. Table 1 and Figure 3 report the result of Model (3).",
   caption = "",
   caption.above = TRUE,
   label = "tab:ind-reg-gov",
   float.pos = "t",
#   file = "figs/table_ind_reg_gov.tex",
   dcolumn = FALSE)
```

Parliament: Table A2
```{r, results = 'asis'}
texreg(list(fits_pmt$m7,
            fits_pmt$m8,
            fits_pmt$m9,
            fits_pmt$m10,
            fits_pmt$m11,
            fits_pmt$m12),
   stars = numeric(0),
   override.se = list( # p-values instead of se
     fits_pmt$m7$p.value,
     fits_pmt$m8$p.value,
     fits_pmt$m9$p.value,
     fits_pmt$m10$p.value,
     fits_pmt$m11$p.value,
     fits_pmt$m12$p.value),
   include.ci = FALSE,
   include.rsquared = FALSE,
   include.adjrs = TRUE,
   include.nobs = TRUE,
   include.fstatistic = FALSE,
   include.rmse = FALSE,
   digits = 3,
   custom.model.names = paste("Model", 7:12),
   custom.note = "Note: p-values calculated with robust standard errors clustered by country are in parentheses. p-values refer to two-tailed tests. All models include country and year fixed-effects. Table 1 and Figure 3 report the result of Model (9).", 
   caption = "",
   caption.above = TRUE,
   label = "tab:ind-reg-pmt",
   float.pos = "t",
#   file = "figs/table_ind_reg_pmt.tex",
   dcolumn = FALSE)
```

EMB: Table A3
```{r, results = 'asis'}
texreg(list(fits_emb$m13,
            fits_emb$m14,
            fits_emb$m15,
            fits_emb$m16,
            fits_emb$m17,
            fits_emb$m18),
   stars = numeric(0),
   override.se = list( # p-values instead of se
     fits_emb$m13$p.value,
     fits_emb$m14$p.value,
     fits_emb$m15$p.value,
     fits_emb$m16$p.value,
     fits_emb$m17$p.value,
     fits_emb$m18$p.value),
   include.ci = FALSE,
   include.rsquared = FALSE,
   include.adjrs = TRUE,
   include.nobs = TRUE,
   include.fstatistic = FALSE,
   include.rmse = FALSE,
   digits = 3,
   custom.model.names = paste("Model", 13:18),
   custom.note = "Note: p-values calculated with robust standard errors clustered by country are in parentheses. p-values refer to two-tailed tests. All models include year fixed-effects. Table 1 and Figure 3 report the result of Model (15).",
   caption = "",
   caption.above = TRUE,
   label = "tab:ind-reg-emb",
   float.pos = "t",
#   file = "figs/table_ind_reg_emb.tex",
   dcolumn = FALSE)
```


# Regression: Effect of blatant fraud

## Models
```{r}
f_gov <- alist(
  m19 = GOV ~ fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud +
    age + female + unemployment + urban,
  m20 = GOV ~ fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud +
    age + female + unemployment + urban) |> 
  map(.f = formula)

f_gov_f <- alist(
  m19 = GOV ~ 0 + fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud +
    age + female + unemployment + urban +
    cowcode + year,
  m20 = GOV ~ 0 + fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud +
    age + female + unemployment + urban +
    cowcode + year) |> 
  map(.f = formula)

f_pmt <- alist(
  m21 = PMT ~ fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud +
    age + female + unemployment + urban,
  m22 = PMT ~ fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud +
    age + female + unemployment + urban) |> 
  map(.f = formula)


f_pmt_f <- alist(
  m21 = PMT ~ 0 + fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud +
    age + female + unemployment + urban +
    cowcode + year,
  m22 = PMT ~ 0 + fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud +
    age + female + unemployment + urban +
    cowcode + year) |> 
  map(.f = formula)

f_emb <-  alist(
  m23 = ETN ~ fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud +
    age + female + unemployment + urban,
  m24 = ETN ~ fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud +
    age + female + unemployment + urban) |> 
  map(.f = formula)

f_emb_f <-  alist(
  m23 = ETN ~ 0 + fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_polyarchy_excfraud +
    age + female + unemployment + urban +
    year,
  m24 = ETN ~ 0 + fraud2 * dif_days10 + 
    nelda17 + nelda18 + 
    gdpgrowth_new + gdppc_change + margin_pty_lag + 
    v2x_libdem_excfraud +
    age + female + unemployment + urban +
    year) |> 
  map(.f = formula)
```

Outcome: Government trust
```{r}
fits_fgov_f <- f_gov_f |> 
  map(.f = lm,
      data = MYDout2)

fits_fgov <- f_gov |> 
  map(.f = lm_robust, 
      data = MYDout2, 
      clusters = cowcode,
      fixed_effects = ~ cowcode + year,
      se_type = "stata")
```


Outcome: Parliament trust
```{r}
fits_fpmt_f <- f_pmt_f |> 
  map(.f = lm,
      data = MYDout2)

fits_fpmt <- f_pmt |> 
  map(.f = lm_robust, 
      data = MYDout2, 
      clusters = cowcode,
      fixed_effects = ~ cowcode + year,
      se_type = "stata")
```

Outcome: EMB trust
```{r}
fits_femb_f <- f_emb |> 
  map(.f = lm,
      data = MYDout2)

fits_femb <- f_emb |> 
  map(.f = lm_robust, 
      data = MYDout2, 
      clusters = cowcode,
      fixed_effects = ~ year,
      se_type = "stata")
```

## Interaction effect

Figure 4: The marginal effect of electoral fraud on popular trust in the government (left panel), parliament (middle), and EMBs (right), conditional on the number of days elapsed since the last election.
```{r}
intf_gov1 <- meplot(fits_fgov_f$m19, 
   vcov = vcov(fits_fgov$m19), 
   data = MYD,
   var1 = "fraud2", 
   var2 = "dif_days10",
   int = "fraud2:dif_days10") +
  labs(x = "Days since last election",
       y = "Effect of electoral fraud",
       title = "Government") +
  ylim(-3, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intf_gov1)

intf_pmt1 <- meplot(fits_fpmt_f$m21, 
   vcov = vcov(fits_fpmt$m21), 
   data = MYD,
   var1 = "fraud2", 
   var2 = "dif_days10",
   int = "fraud2:dif_days10") +
  labs(x = "Days since last election",
       y = "",
       title = "Parliament")  +
  ylim(-3, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intf_pmt1)

intf_emb1 <- meplot(fits_femb_f$m23, 
   vcov = vcov(fits_femb$m23), 
   data = MYD,
   var1 = "fraud2", 
   var2 = "dif_days10",
   int = "fraud2:dif_days10") +
  labs(x = "Days since last election",
       y = "",
       title = "EMB") +
  ylim(-3, 1.25) +
  scale_x_continuous(breaks = seq(0, 200, by = 50),
                     labels = seq(0, 200, by = 50) * 10)
#plot(intf_emb1)

print(intf_gov1 | intf_pmt1 | intf_emb1)
```

Regression results: Table 2
```{r, results = 'asis'}
texreg(list(fits_fgov$m19,
            fits_fpmt$m21,
            fits_femb$m23),
   stars = numeric(0),
   override.se = list( # p-values instead of se
     fits_fgov$m19$p.value,
     fits_fpmt$m21$p.value,
     fits_femb$m23$p.value
   ),
   include.ci = FALSE,
   include.rsquared = FALSE,
   include.adjrs = TRUE,
   include.nobs = TRUE,
   include.fstatistic = FALSE,
   include.rmse = FALSE,
   digits = 3,
   custom.model.names = c("gov", "pmt", "emb"),
   custom.note = "Note: p-values calculated with robust standard errors clustered by country are in parentheses. p-values refer to two-tailed tests. The first two models include year-fixed effects and all models include country-fixed effects.",
   caption = "",
   caption.above = TRUE,
   label = "tab:ols-res-main",
   float.pos = "t",
#   file = "figs/table_ols_res_fraud_main.tex",
   dcolumn = FALSE)
```

Regression results: Table A4
```{r, results = 'asis'}
texreg(list(fits_fgov$m19,
            fits_fgov$m20,
            fits_fpmt$m21,
            fits_fpmt$m22,
            fits_femb$m23,
            fits_femb$m24),
   stars = numeric(0),
   override.se = list( # p-values instead of se
     fits_fgov$m19$p.value,
     fits_fgov$m20$p.value,
     fits_fpmt$m21$p.value,
     fits_fpmt$m22$p.value,
     fits_femb$m23$p.value,
     fits_femb$m24$p.value
   ),
   include.ci = FALSE,
   include.rsquared = FALSE,
   include.adjrs = TRUE,
   include.nobs = TRUE,
   include.fstatistic = FALSE,
   include.rmse = FALSE,
   digits = 3,
   custom.model.names = paste("Model", 19:24),
   custom.note = "Note: p-values calculated with robust standard errors clustered by country are in parentheses. p-values refer to two-tailed tests. Models for government trust and parliament trust include country and year fixed-effects. The rest includes the year fixed-effect. Table 2 and Figure 4 report the results of Models (19), (21), and (23).",
   caption = "",
   caption.above = TRUE,
   label = "tab:ols-res-fraud",
   float.pos = "t",
#   file = "figs/table_ind_reg_fraud.tex",
   dcolumn = FALSE)
```

We conducted multilevel regression (results shown in Table A5) in Stata.


